!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


module mod_octave_io
!General comments: write variables in octave text form. It is assumed
! that the file is already open and connected to unit fu.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_octave_io_constructor, &
   mod_octave_io_destructor,  &
   mod_octave_io_initialized, &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   real_format,    &
   locate_var

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! private members
 character(len=*), parameter :: &
   real_format = 'e24.15e3' ! format: _+0.***************E+***
!   real_format = 'e42.33e3' ! format: _+0.** 33 **E+*** : real(16)

! Module variables

 ! public members
 logical, protected ::               &
   mod_octave_io_initialized = .false.

 interface write_octave
   module procedure write_scal_r, write_vect_r, write_mat_r,  &
                    write_3dim_r, write_4dim_r, write_5dim_r, &
                    write_scal_i, write_vect_i, write_mat_i,  &
                    write_3dim_i, write_4dim_i, write_5dim_i, &
                    write_scal_l, write_vect_l, write_mat_l,  &
                    write_3dim_l, write_4dim_l, write_5dim_l, &
                    write_scal_c, write_mult_c
 end interface

 interface read_octave
   module procedure read_scal_r, read_vect_r, read_mat_r,  &
                    read_3dim_r, read_4dim_r, read_5dim_r, &
                    read_scal_i, read_vect_i, read_mat_i,  &
                    read_3dim_i, read_4dim_i, read_5dim_i, &
                    read_scal_l, read_vect_l, read_mat_l,  &
                    read_3dim_l, read_4dim_l, read_5dim_l
 end interface

 interface read_octave_al
   module procedure read_vect_al_r, read_mat_al_r,  &
!                    read_3dim_r, read_4dim_r, read_5dim_r, &
                    read_vect_al_i, read_mat_al_i!,  &
!                    read_3dim_i, read_4dim_i, read_5dim_i, &
!                    read_scal_l, read_vect_l, read_mat_l,  &
!                    read_3dim_l, read_4dim_l, read_5dim_l
 end interface

 character(len=*), parameter :: &
   this_mod_name = 'mod_octave_io'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_octave_io_constructor()

  character(len=100) message
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_octave_io_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   write(message,'(A,A)') 'Octave output format: ', real_format
   call info(this_sub_name,this_mod_name,message)

   mod_octave_io_initialized = .true.
 end subroutine mod_octave_io_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_octave_io_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_octave_io_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_octave_io_initialized = .false.
 end subroutine mod_octave_io_destructor

!-----------------------------------------------------------------------
 
 subroutine write_scal_r(x,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: x
  character(len=*), intent(in) :: var_name
 
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

  write(fu,'(a,a)')    '# name: ',var_name
  write(fu,'(a)')      '# type: scalar'

  write(out_format,'(a,a,a)') '(',real_format,')'
  write(fu,out_format) x
  
 end subroutine write_scal_r
 
!-----------------------------------------------------------------------

 subroutine read_scal_r(x,var_name,fu,norewind)
  integer, intent(in) :: fu
  character(len=*), intent(in) :: var_name
  real(wp), intent(out) :: x
  logical, intent(in), optional :: norewind
 
  integer :: ierr
  character(len=1000) :: text
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

  ! the following subroutine positions the file at the '# name' line
  call locate_var(fu,var_name,ierr,norewind)
  if(ierr.eq.0) then
    read(fu,*) text
    read(fu,*) x
  else
    write(text,'(A,A,A)') 'Problem reading variable "',var_name,'".'
    call warning(this_sub_name,this_mod_name,text)
  endif
  
 end subroutine read_scal_r
 
!-----------------------------------------------------------------------
 
 subroutine write_vect_r(V,rc,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: V(:)
  character(len=*), intent(in) :: rc, var_name

  integer :: nrow, ncol
  character(len=17) :: out_format
  character(len=100) :: message(2)
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   if(rc.eq.'c') then
     nrow = product(shape(V))
     ncol = 1
   elseif(rc.eq.'r') then
     nrow = 1
     ncol = product(shape(V))
   else
     write(message(1),'(A)') &
       'Variable rc must be either "r" or "c".'
     write(message(2),'(A,A,A)') &
       'Present value rc="',rc,'".'
     call warning(this_sub_name,this_mod_name,message)
     return
   endif
     
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   ! Note: the repeat count must be strictly positive; on the other
   ! hand it doesn't have to match the data size, so if it's zero we
   ! should set it at least to one
   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, real_format, ')'
   write(fu,out_format) V
  
 end subroutine write_vect_r
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_r(V,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, nrow, ncol
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'
 
   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,'(A)') message(1) ! exclude "# type: matrix"
     read(fu,'(A7,I10)') message(1), nrow ! "# rows:"
     read(fu,'(A10,I10)') message(1), ncol ! "# columns:"
     if((nrow.ne.1).and.(ncol.ne.1)) then
       write(message(1),'(A)') &
         'Variable "',var_name,'" is not a vector.'
       write(message(2),'(A,I10,A,I10)') &
         '  dimension: nrow=',nrow,', ncol=',ncol
       call warning(this_sub_name,this_mod_name,message(1:2))
     elseif(product(shape(V)).ne.(nrow*ncol)) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10)') &
         'dimension in file: ',ncol,nrow
       write(message(3),'(A,I10)') &
         'internal dimension: ',shape(V)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) V
     endif
   endif
  
 end subroutine read_vect_r
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_al_r(V,var_name,fu,norewind)
  integer, intent(in) :: fu
  real(wp), allocatable, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, inrow, incol
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"
     ! A matrix with only one element is saved by octave as a scalar,
     ! which requires a special case.
     if(trim(adjustl(type)).eq."scalar") then
       inrow = 1
       incol = 1
     else
       read(fu, '(A7,I10)') message, inrow ! "# rows:"
       read(fu,'(A10,I10)') message, incol ! "# columns:"
     endif
     allocate( V(inrow*incol) )
     read(fu,*) V
   endif
  
 end subroutine read_vect_al_r
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_r(M,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: nrow, ncol
  character(len=17) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   nrow = size(M,1)
   ncol = size(M,2)
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   ! Note: the repeat count must be strictly positive; on the other
   ! hand it doesn't have to match the data size, so if it's zero we
   ! should set it at least to one
   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, real_format, ')'
   if(ncol.ne.0) & ! otherwise nothing to write
     write(fu,out_format) transpose(M)
  
 end subroutine write_mat_r
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_r(M,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, nrow, ncol, i
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'
 
   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,'(A)') message(1) ! exclude "# type: matrix"
     read(fu,'(A7,I10)') message(1), nrow ! "# rows:"
     read(fu,'(A10,I10)') message(1), ncol ! "# columns:"
     if( (size(M,1).ne.nrow).or.(size(M,2).ne.ncol) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10)') &
         'dimension in file: ',ncol,nrow
       write(message(3),'(A,I10,I10)') &
         'internal dimension: ',shape(M)
       call warning(this_sub_name,this_mod_name,message)
     else
       do i=1,nrow
         read(fu,*) M(i,:)
       enddo
     endif
   endif
  
 end subroutine read_mat_r
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_al_r(M,var_name,fu,norewind)
  integer, intent(in) :: fu
  real(wp), allocatable, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, inrow, incol, i
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"
     ! A matrix with only one element is saved by octave as a scalar,
     ! which requires a special case.
     if(trim(adjustl(type)).eq."scalar") then
       inrow = 1
       incol = 1
     else
       read(fu, '(A7,I10)') message, inrow ! "# rows:"
       read(fu,'(A10,I10)') message, incol ! "# columns:"
     endif
     allocate( M(inrow,incol) )
     do i=1,inrow
       read(fu,*) M(i,:)
     enddo
   endif
  
 end subroutine read_mat_al_r
 
!-----------------------------------------------------------------------
 
! Unfortunately, an assumed size array T(*), which would cover all the
! 3D, 4D, xD cases, would not match the generic interface. Thus, it
! would be necessary to call directly the write_multidim subroutine.
! The following is nicer, provided we don't want really large rank
! arrays.
 subroutine write_3dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(3)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',3
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_3dim_r
 
!-----------------------------------------------------------------------

 subroutine read_3dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(3)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3)
       write(message(3),'(A,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_3dim_r
 
!-----------------------------------------------------------------------

 subroutine write_4dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(4)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',4
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_4dim_r
 
!-----------------------------------------------------------------------

 subroutine read_4dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(4)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)).or. &
         (size(T,4).ne.dims(4)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3),dims(4)
       write(message(3),'(A,I10,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_4dim_r
 
!-----------------------------------------------------------------------

 subroutine write_5dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(5)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',5
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_5dim_r
 
!-----------------------------------------------------------------------

 subroutine read_5dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(5)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)).or. &
         (size(T,4).ne.dims(4)).or. &
         (size(T,5).ne.dims(5)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3), &
                               dims(4),dims(5)
       write(message(3),'(A,I10,I10,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_5dim_r

!-----------------------------------------------------------------------

 subroutine write_scal_i(x,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: x
  character(len=*), intent(in) :: var_name
 
   call write_scal_r(real(x,wp),var_name,fu)

 end subroutine write_scal_i

!-----------------------------------------------------------------------

 subroutine read_scal_i(x,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, intent(out) :: x
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind
 
  real(wp) :: x_r

   call read_scal_r(x_r,var_name,fu,norewind)
   x = int(x_r)

 end subroutine read_scal_i

!-----------------------------------------------------------------------
 
 subroutine write_vect_i(V,rc,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: V(:)
  character(len=*), intent(in) :: rc, var_name

   call write_vect_r(real(V,wp),rc,var_name,fu)
  
 end subroutine write_vect_i
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_i(V,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  real(wp) :: V_r(size(V,1))

   call read_vect_r(V_r,var_name,fu)
   V = int(V_r)
  
 end subroutine read_vect_i
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_al_i(V,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, allocatable, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  real(wp), allocatable :: V_r(:)

   call read_vect_al_r(V_r,var_name,fu,norewind)
   allocate( V(size(V_r)) )
   V = int(V_r)
  
 end subroutine read_vect_al_i
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_i(M,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

   call write_mat_r(real(M,wp),var_name,fu)

 end subroutine write_mat_i
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_i(M,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: M_r(size(M,1),size(M,2))

   call read_mat_r(M_r,var_name,fu)
   M = int(M_r)
  
 end subroutine read_mat_i
 
!-----------------------------------------------------------------------

 subroutine read_mat_al_i(M,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, allocatable, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  real(wp), allocatable :: M_r(:,:)

   call read_mat_al_r(M_r,var_name,fu,norewind)
   allocate( M(size(M_r,1),size(M_r,2)) )
   M = int(M_r)
  
 end subroutine read_mat_al_i
 
!-----------------------------------------------------------------------

 subroutine write_3dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

   call write_3dim_r(real(T,wp),var_name,fu)

 end subroutine write_3dim_i
 
!-----------------------------------------------------------------------

 subroutine read_3dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3))

   call read_3dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_3dim_i
 
!-----------------------------------------------------------------------

 subroutine write_4dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

   call write_4dim_r(real(T,wp),var_name,fu)

 end subroutine write_4dim_i
 
!-----------------------------------------------------------------------

 subroutine read_4dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3),size(T,4))

   call read_4dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_4dim_i
 
!-----------------------------------------------------------------------

 subroutine write_5dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

   call write_5dim_r(real(T,wp),var_name,fu)

 end subroutine write_5dim_i
 
!-----------------------------------------------------------------------

 subroutine read_5dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   call read_5dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_5dim_i
 
!-----------------------------------------------------------------------

 subroutine write_scal_l(x,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: x
  character(len=*), intent(in) :: var_name
 
  integer :: x_i

   if(x.eqv..true.) then
     x_i = 1
   else
     x_i = 0
   endif
   call write_scal_i(x_i,var_name,fu)

 end subroutine write_scal_l

!-----------------------------------------------------------------------

 subroutine read_scal_l(x,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: x
  character(len=*), intent(in) :: var_name
 
  integer :: x_i

   call read_scal_i(x_i,var_name,fu)
   if(x_i.eq.1) then
     x = .true.
   else
     x = .false.
   endif

 end subroutine read_scal_l

!-----------------------------------------------------------------------
 
 subroutine write_vect_l(V,rc,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: V(:)
  character(len=*), intent(in) :: rc,var_name

  integer :: V_i(size(V,1))

   where(V)
     V_i = 1
   elsewhere
     V_i = 0
   endwhere
   call write_vect_i(V_i,rc,var_name,fu)
  
 end subroutine write_vect_l
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_l(V,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  integer :: V_i(size(V,1))

   call read_vect_i(V_i,var_name,fu)
   where(V_i.eq.1)
     V = .true.
   elsewhere
     V = .false.
   endwhere
  
 end subroutine read_vect_l
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_l(M,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: M_i(size(M,1),size(M,2))

   where(M)
     M_i = 1
   elsewhere
     M_i = 0
   endwhere
   call write_mat_i(M_i,var_name,fu)
  
 end subroutine write_mat_l
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_l(M,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: M_i(size(M,1),size(M,2))

   call read_mat_i(M_i,var_name,fu)
   where(M_i.eq.1)
     M = .true.
   elsewhere
     M = .false.
   endwhere
  
 end subroutine read_mat_l
 
!-----------------------------------------------------------------------
 
 subroutine write_3dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_3dim_i(T_i,var_name,fu)
  
 end subroutine write_3dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_3dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3))

   call read_3dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_3dim_l
 
!-----------------------------------------------------------------------
 
 subroutine write_4dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_4dim_i(T_i,var_name,fu)
  
 end subroutine write_4dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_4dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4))

   call read_4dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_4dim_l
 
!-----------------------------------------------------------------------
 
 subroutine write_5dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_5dim_i(T_i,var_name,fu)
  
 end subroutine write_5dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_5dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   call read_5dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_5dim_l
 
!-----------------------------------------------------------------------

 subroutine write_scal_c(s,var_name,fu)
  integer, intent(in) :: fu
  character(len=*), intent(in) :: s
  character(len=*), intent(in) :: var_name
 
  write(fu,'(a,a)')    '# name: ',var_name
  write(fu,'(a)')      '# type: sq_string'
  write(fu,'(a)')      '# elements: 1'
  write(fu,'(a,i7)')   '# length: ',len(trim(s))
  write(fu,'(a)')      trim(s)

 end subroutine write_scal_c

!-----------------------------------------------------------------------

 subroutine write_mult_c(s,var_name,fu)
 ! Write a character array on multiple lines; each line is trimmed.
  integer, intent(in) :: fu
  character(len=*), intent(in) :: s(:)
  character(len=*), intent(in) :: var_name
 
  integer :: lt, i

  write(fu,'(a,a)')    '# name: ',var_name
  write(fu,'(a)')      '# type: sq_string'
  write(fu,'(a)')      '# elements: 1'
  ! determine the total length: first count the newline characters
  lt = max(size(s)-1,0)
  do i=1,size(s)
    lt = lt + len(trim(s(i)))
  enddo
  write(fu,'(a,i7)')   '# length: ',lt
  do i=1,size(s)
    write(fu,'(a)')    trim(s(i))
  enddo

 end subroutine write_mult_c

!-----------------------------------------------------------------------
 
 subroutine locate_var(fu,var_name,ierr,norewind)
 ! Position the file at the first occurrence of string
 ! Use  norewind.eqv..true.  to avoid rewind
  integer, intent(in) :: fu
  character(len=*), intent(in) :: var_name
  integer, intent(out) :: ierr
  logical, intent(in), optional :: norewind

  logical :: found, rew
  character(len=1000) :: record
  character(len=len(var_name)+8) :: record_var_name
  character(len=*), parameter :: &
    this_sub_name = 'locate_var'

   found = .false.
   record_var_name = '# name: ' // var_name! // '\n'

   ! Check whether rewind is required
   rew = .true. ! default is to rewind the file
   if(present(norewind)) rew = .not.norewind
   if(rew) rewind(fu)

   do
    if(found) exit
     read(fu,'(a)',iostat=ierr) record
     if(ierr.ne.0) then ! either error or end of file
       found = .true.
     elseif(record.eq.record_var_name) then
       found = .true.
       ierr = 0
     endif
   enddo

 end subroutine locate_var

!-----------------------------------------------------------------------

end module mod_octave_io

